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ABSTRACT 

We present virial models for the global evolution of giant molecular clouds. Focusing on the presence 
of an accretion flow, and accounting for the amount of mass, momentum, and energy supplied by 
accretion and star formation feedback, we are able to follow the growth, evolution, and dispersal of 
individual giant molecular clouds. Our model clouds reproduce the scaling relations observed in both 
galactic and extragalactic clouds. We find that accretion and star formation contribute roughly equal 
amounts of turbulent kinetic energy over the lifetime of the cloud. Clouds attain virial equilibrium and 
grow in such a way as to maintain roughly constant surface densities, with typical surface densities of 
order 50 - 200 M pc -2 , in good agreement with observations of giant molecular clouds in the Milky 
Way and nearby external galaxies. We find that as clouds grow, their velocity dispersion and radius 
must also increase, implying that the linewidth-size relation constitutes an age sequence. Lastly, we 
compare our models to observations of giant molecular clouds and associated young star clusters in 
the LMC and find good agreement between our model clouds and the observed relationship between 
H II regions, young star clusters, and giant molecular clouds. 

Subject headings: stars: formation — ISM: clouds — ISM: evolution — turbulence 



1. INTRODUCTION 

Giant molecular clouds (GMCs) ar e the primary reser- 
voir o f molecular g as in t he galaxy (iWilliams fe McKed 
[19971: IRosolowskvl [20051: IStark fc Led 12006ft . Since 
the surface density of star formation shows a strong 
correlation with the surface density of molecular gas 
(iWong fc Blitdl200a iKennicutt et"aJ]|2007t iBigiel et all 
I2008t iSchruba et al J 1201 If) . GMCs rnust also be the pri- 
mary site of star formation in the Milky Way. However, 
recent high-resolution observations have shown that the 
Kennicutt-Schmidt law breaks down when the resolution 
of an o bservation is finer than the typical length s cales of 
GMCs (IQnodera et al.ll20Tot ISchruba et alJl20Toh . Thus, 
in order to develop a detailed theoretical understanding 
of the relationship between star formation and molecu- 
lar gas, it is necessary to first understand the formation, 
evolution, and destruction of giant molecular clouds. 

One stumbling block in this effort is the substan- 
tial disagreement in the literature regarding both the 
formatio n mechanism and typica l lifetimes of GMCs 
see e.g. iGoldreich fc Kwanl 119741: iZuckerman fc Evansl 



1998 : iStone et al 



Elmegreen fc Scalo 



197llBlitz fc Shul)1980t iBallesteros-Paredes et all 119991 : 
McKee fc Ostrikerl l2007t iMurravl I2011L and references 



therein). Some authors suggest that GMCs form out 
of bound atomic gas as a result of gravitational instabil- 
ity (IKim et al.l2lMl200ltlKirn fc Ostrikerl2006HLi et all 
I2006at iTasker fc Tanl 120031 surviving as rou ghly virial- 
ized o bjects for many cloud dynamical times (jTan et all 
120061: iTamburro et all 120081 ). In support of this pic- 
ture is the observation that massive clouds are found 
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to be marginally bou nd, with typical v i rial param- 
eters of order unity dHeyer et all 120011 : IRosolowskvl 
[20071 iRoman-Duval et all l20io[h bmce supersonic 
isothermal turbulence is found to dec ay via radiative 
shock s in one or t wo cr os sing times dMac Low et all 

iMac Low fc Klessenl 120041 : 
this model must invoke a 
mechanism to drive supersonic motions for the life- 
time of a cloud, which could be several crossing times. 
Possible turbulen t driving mechanisms include proto- 
stellar outflows (iNorman fc SilkJ [19801 : IMcKed 119891: 
iLi fc Nakamural 120061 iLi et all 12010ft H II regions 
(jMatznerl l2002t iKrumholz fc Matznerl 120091 h ereafter 
KM09) , supernovae (|Mac Low fc Klessenl I200l. or, as 
investigated here, mass accretion ((Klessen fc Hennebelld 
l2010tlVazquez-Semadeni et al.ll2010l) . 

Accretion driven turbulence in molecular clouds has 
received little attention in t he literature. However, as 
IKlessen fc Hennebelld (|2010h point out, the kinetic en- 
ergy of accreted material can power the turbulent mo- 
tions observed in molecular clouds with energy conver- 
sion efficiencies of only a few percent. While there has 
been no systematic study of the kinetic energy budget 
of a molecular cloud formed via gravitational instability, 
this problem has been examined in the context of the for- 
mation of pro t ogalax ies at high rcdshift. In one example, 
iWise fc Abell (|2007t) analyzed simulations of virializing 
high redshift minihalos, tracking the thermal, kinetic, 
and gravitational potential energy of gas in protogalac- 
tic dark matter halos. In their models, which included 
a nonequilibrium cooling model, gas collapsed onto the 
halo and cooled quickly, causing turbule nt velocities to 
becom e supersonic. As pointed out by IWang fc Abell 
(|200l . this means that the virialization process is a lo- 
cal one: gravitational potential energy can be converted 
directly into supersonically turbulent motions character- 
ized by a volume-filling network of shocks. The turbu- 
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lence in turn provides much of the kinetic support for the 
newly virializcd gaseous component of the dark matter 
halo. If a similar mechanism is at work as gas cools and 
collapses onto GMCs then gravitational potential energy 
alone may be sufficient to power turbulence in GMC. 

The most detailed simulations of simultaneously ac- 
creting and star form ing giant molecular clou d s wer e 
recently completed by IVazquez-Semadeni et al.l (|2010l ). 
These numerical models included a simplified subgrid 
prescription for stellar feedback by the ionizing radia- 
tion of newborn star clusters and focused on the bal- 
ance between accretion and feedback in clouds formed 
via thermal instability in colliding flows. Throughout the 
course of the simulation, dense molecular gas condensed 
out of a warm atomic envelope, allowing a study of the 
interplay between accretion and feedback in the simu- 
lated clouds. The resulting clouds were able to attain a 
state of quasi-virial equilibrium, in which the supply of 
gas from the ambient medium balanced the formation of 
stars and ejection of gas from H n regions. Due to the 
idealized nature of the subgrid star formation feedback 
prescription, in which all H n regions were powered by a 
cluster with the same ionizing luminosity, star formation 
feedback was unable to act on the cloud as a whole but 
could reduce the global star formation rate by destroy- 
ing overdensities. Since the simulation did not include 
star clusters with large ionizing luminosities, the cloud 
as a whole could not be destroyed and star formation 
would have eventually consumed all of the gas had the 
simulation not been cut off. Even though the simulations 
employed a highly idealized star formation prescription, 
the computations still required substantial resources to 
complete and only allowed insight into the evolution of 
a single cloud. It seems that a computationally inexpen- 
sive model that includes a somewhat more sophisticated 
treatment of star formation feedback is called for. 

In this work, we model the global evolution of giant 
molecular clouds from their birth as low-mass seed clouds 
to their dispersal after a phase of massive star formation. 
This is done using an updated version of the semia nalyt- 
ical model of lKrumholz. Matzner k McKee I (|2006l here- 
after Paper I). Using a virial formalism, we compute the 
global dynamical evolution of a single cloud while simul- 
taneously tracking its energy budget. Model clouds form 
stars, launch H n regions and undergo accretion from 
their environments. With these models, we are able to 
investigate the role accretion plays in maintaining turbu- 
lence in molecular clouds and directly compare to obser- 
vations of GMCs in the Milky Way and nearby external 
ga laxies. This work is comp l ement ary to the simulations 
of IVazquez -Semad eni et all (|2010f ) . since our simplified 
global models allow us to survey a large variety of GMCs 
at little computational cost while including a much more 
sophisticated star formation feedback prescription. We 
are able to capture model clouds with masses compara- 
ble to the most massive clouds observed in the Milky 
Way and nearby galaxies, allowing us to simulate the 
sites of the majority of star formation in these syst ems 
([Williams k McKeejfl997l iFukui k Kawamurall2010l ). 

We proceed by describing the formulation and imple- 
mentation of our GMC model in §[5] Next, in §[31 we 
test our implementation of accretion. Following this, in 
§[4] we perform full simulations and describe the gneral 
features of our simulated clouds. In §[5]we make compar- 



isons to observations, focusing on the scaling relations 
observed to hold for GMCs as well as the high quality 
multiwavclcngth observations available for GMCs in the 
LMC. Lastly, in §[6l we discuss the limitations inherent in 
the simplifying assumptions we make to derive the cloud 
evolution equations. 

2. GOVERNING EQUATIONS 

The GMC evolution model described below allows us 
to solve for the time-evolution of the global properties of 
model molecular clouds. In contrast with previous work, 
we follow the flow of gas as it condenses out of the dif- 
fuse gas in the envelope surrounding the GMC and falls 
onto the cloud. Employing simplifying assumptions as 
well as the results of simulations of compressible MHD 
turbulence, we derive a set of coupled ordinary differ- 
ential equations that govern the time evolution of the 
cloud's mass, radius, and velocity dispersion. Combin- 
ing the governing evolution equations with a set of ini- 
tial conditions, model parameters, and a model for the 
time dependence of the mass accretion rate based on the 
gravitational collapse of the GMC envelope, we can solve 
for the time evolution of the cloud. Below, we give an 
overview of the model, discuss the formulation of our nu- 
merical scheme, describe our parameter choices, and give 
a brief description of our treatment of star formation and 
our model for the GMCs gas supply. 

2.1. Model Overview 

The model we employ below is based on the global 
GMC model of Paper I, itself a generalization of 
the the globa l model for low mass star formation of 
McKed (119891) . the Eul erian virial theorem (EVT) of 
McKee fe Zweibel (119921) . and the model for star forming 
clumps of lMatznerl (|200lh . Employing a virial formalism, 
we account for the dynamics and energy budget of gas 
contained within an Eulerian volume, V v i r . We separate 
the gas within V v i r into three species: virial material, a 
gaseous reservoir, and a photoionized wind. A schematic 
representation of the components of our model is pre- 
sented in Figured] 

By design, each of the three components has a straight- 
forward physical interpretation. The first component, 
which we label virial material, consists of two physi- 
cally distinct subcomponents: a molecular cloud and 
a warm atomic envelope that encloses the cloud. The 
cloud is assumed to be cold (~ 10 K), molecular, and 
contained within a spherical volume of radius R c \. The 
ambient medium is composed of warm (~ 10 3 K) and 
diffuse atomic gas that encloses the cloud and extends 
beyond the virial volume. The second component is a 
gaseous reservoir, which we assume is composed of cold 
(~ 10 2 K) neutral material that flows onto the cloud at 
free-fall from beyond the virial volume. The last compo- 
nent is an ionized wind made up of hot (~ 10 4 K) ionized 
gas ejected from the ionization fronts of blister-type H n 
regions. All three components are allowed to mutually 
interpenetrate. We restrict interaction between the com- 
ponents to the transfer of mass between the accretion 
flow and cloud as well as between the cloud and wind. 
Since the envelope and cloud are not allowed to inter- 
penetrate, we formally group the envelope and cloud to- 
gether; this somewhat artificial choice significantly sim- 
plifies the virial analysis. 
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Figure 1. A schematic overview of the GMC model. A molec- 
ular cloud (black) is embedded in a warm atomic envelope (dark 
blue). Cool atomic gas (light blue) flows onto the cloud, where it 
condenses, recombines into molecules, and mixes with the cloud. 
Newborn OB associations (blue stars) drive H II regions (orange) 
and eject winds ionized winds back into the ambient medium. 

Wc make use of two simplifying assumptions regarding 
the distribution and flow of the virial material. First, 
we assume that the virial material follows a spherically 
symmetric, smoothly varying density profile. Second, we 
assume that the the cloud is homologous: the cloud ex- 
pands, contracts, accretes, and sheds mass in such a way 
as to always maintain the same smooth density profile. 
We assume a density profile of the form 



P( r ) =Po[-5- 



r 

Rci 



for r < R c 



(1) 



where po is the density at the edge of the cloud and k p 
is assumed to be unity. This choice is consiste nt with 
the Larson scaling relations observed in galactic (|Larsonl 
fl98lt [Solornon et all 119871: iHever et al.ll2001l 120091) and 
extragalacti c (IMizuno et al.l 2001 [ Engargiola et alJ2003t 
Rosolowskvl 120071: IBolatto et al.l 120081: iHughes et al l 
201C0 GMCs. Our assumed cloud density profile ef- 



fectively averages over clumpy and filamentary inter- 
nal structure and oblate shapes of observed clouds. At 
r = R c i, the density of the ambient medium is assumed 
to smoothly transition from pc to /3 a mb in a thin bound- 
ary layer. We assume that the density of the ambient 
medium is negligible compared to the density of gas 
in the cloud, p a mb "C po- The density of the ambient 
medium remains p am b out to the virial radius. 

Beyond assuming a density profile, we must also spec- 
ify the velocity structure of all three components of the 
model. We follow Paper I in assuming that the velocity of 
the virial material can be decomposed into a systematic 
and turbulent component, 



Rc 



Vturb- 



(2) 



We assume that v tur b is randomly oriented with respect 
to position so that turbulent motions carry no net flux 



of matter. We make a similar assumption regarding the 
velocity structure of the reservoir, 

Vres ^res,sys^ Vres,turb* (3) 

The systematic component of v res is due to the gravita- 
tional attraction of material within the virial volume 



g • dr, 



(4) 



while 

r-l 



the random component is such that 
(^cT7y cl PvL,turb^) 1/2 = V3cw Here g is the 
gravitational acceleration and <7 res is the velocity disper- 
sion of gas in the reservoir feeding the accretion flow. 
Since the amount of material in the accretion flow is 
determined by how fast it can fall into the cloud, wc 
must simultaneously determine both the density profile 
and radial velocity of material in the accretion flow (see 
Appendix [B| . Finally, for the wind material, we follow 
Paper I in assuming 

v„ = v + (5) 

where = 2c n f and c„ is the ionized gas sound speed. 
Wc follow iMcKee fc Williamsl (|1997h in choosing c n = 
9.7 km s" 1 . 

2.2. Momentum Equation 

In Appendix \Xl we derive the EVT for a simultane- 
ously evaporating and accreting cloud, 



i/ CI = 2( r-ro) +B + >v-i| 



(pvr 2 ) • dS 



aiM cl R d R c i + -ajMdR^ + aiM cs R cl R cl (6) 



3 k p 

4 - k, 



R c i(M ei v' ej - £M., 



acc L/ csc I 



Here, a\ is a constant of order unity that depends on 
the distribution of material in the cloud, I c \ is the cloud 
moment of inertia, T is the combined turbulent and 
thermal kinetic energy of the cloud, 7o is the energy 
associated with interstellar pressure at the cloud sur- 
face, B is the net magnetic energy due to the pres- 
ence of the cloud, W is the gravitational term (equal 
to the gravitational binding energy in the a bsence of 
an external potential McKce & Zweibel 1992]), the sur- 
face integral is proportional to the rate of change of 
the moment of inertia inside the bounding viral surface, 
M c \ is the cloud mass, R c \ is the cloud radius, M e j is 
the mass ejection rate, Af aC c is the mass accretion rate, 
v csc = {2G[M cl + M^R^/Ra} 1 ' 2 is the escape ve- 
locity at the edge of the cloud, and £ is a dimcnsionless 
factor we compute via equation (|A8|) that depends on the 
depth of the cloud potential well. The quantity £i> osc is 
the accretion rate weighted average infall velocity. Pre- 
cise definitions for T, To, B, and W are given in Paper 
I. 

The EVT of a cloud without accretion or mass loss 
would only contain the terms up to the surface integral. 
The next three terms account for changes in the cloud 
moment of inertia due to changes in the mass of the 
cloud, while the last term accounts for the rate at which 
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the recoil of inflowing and outflowing mass injects mo- 
mentum into the cloud. Inflows and outflows are treated 
separately because material is ejected at a constant ve- 
locity, but is accreted at a velocity that is a function of 
the depth of the potential well of the cloud. The dimcn- 
sionlcss factor £ appears due to this difference. 

The mass of the cloud can only change via mass accre- 
tion or ejection, 



M d = M, 



M a 



(7) 

We assume that ejection of material can only decrease 
the mass of the cloud and accretion can only increase 
the mass of the cloud. Since stars may not follow the 
homologous density profile we assume, we neglect the 
change in the cloud mass due to star formation. We 
expect the error incurred from this assumption will be 
small, since star s make up a small fraction of the ma ss of 
observed clouds (jEvans et al.ll2009t lLada et aT1l2010ft and 
our assumed star formation law converts a small fraction 
of the cloud's mass into stars per free-fall time. 

We follow Paper I in using the assumption of homol- 
ogy and the results of simulations of MHD turbulence to 
evaluate each term in the EVT in terms of constants and 
the dynamical variables R c \, M c \, and a c \. In the end, 
we obtain a second order nonlinear ordinary differential 
equation in R c \, 



aiR c 




(8) 



Here, 



15 - 5k p 
15- 6fc c 



14- (3 - 2k p 



y{x)dx 



(9) 



y(x) is the ratio of the mass of reservoir material to the 
mass of cloud material contained within a normalized 
radius x = r / R c \ (see Appendix [Bj , and ?7b is the ratio 
of the magnetic critical mass to the cloud mass. 

This equation governs the balance of forces acting on 
the cloud as a whole. Each term corresponds to a sin- 
gle physical mechanism that can alter the radial force 
balance. The first two terms are due to thermal and tur- 
bulent pressure support, respectively. The third is due to 
a combination of gravitational compression and magnetic 
support. The fourth is due to the confining interstellar 
pressure. The fifth comes from the exchange of momen- 
tum between the expanding cloud and infalling accretion 
flow. Finally, the last term is due to a combination of 
the recoil from ejected material and the ram pressure of 
accreting material. Although the two parts of the recoil 
term have opposite signs, M e j and M acc have opposite 

signs as well: M j < and M acc > 0. This implies that 
both the recoil due to launching wind material and the 
ram pressure of accreting reservoir material tend to con- 
fine the cloud. 

Letting M c \$, R c \,a, and (Tc^o be the cloud mass, radius, 
and velocity dispersion at t = and defining the initial 



cloud crossing time, i CI% o = -Rci.o/cci.o, we can define the 
dimensionless variables M = M c \/M c \ y o, R = R c \/R c \ t o, 
= Cci/cci.O: an d t = t/t CI _o. Letting primes denote dif- 
ferentiation with respect to r, we can write equation ([8]) 
in dimensionless form 



R 



where 



„ 3.9cr 2 4- 3M- £ 



ajR 
M'R 



r\G- 



a'M 



R*_ 

R 2 VP M 



M (fMR) 1 / 2 



Mo = <7cl,o/Ccl 



(10) 



(11) 



is the initial turbulent Mach number and we define the 
dimensionless constants 



VG 



VP 



VE 



VA 



3(1-/71) 

a/"vir,0 
47r^ 10 P amb 

a/ilicLoofj 




where 



Ckvir.O 



GA/ Cl . 



(12) 
(13) 
(14) 
(15) 

(16) 



is the initial non thermal virial parameter 
(jBertoldi fc McKed 119921 ). These constants are set 
by the ratios of various forces acting on the initial 
state of the cloud, tjq is proportional to the ratio of 
the initial magnetic forces to the initial gravitational 
force, and rjp, and t\a are the ratios of the ambient 
pressure force, the mass ejection recoil force, and the 
initial accretion ram pressure force to the initial internal 
turbulent forces, respectively. 

Comparing equation (fTU|) with the corresponding equa- 
tion given in Paper I, we see that two new terms pro- 
portional to M acc have appeared. In practice, we find 
that, of the two terms, the one proportional to t\a domi- 
nates, implying that the primary direct impact of accre- 
tion on the radial force balance of the cloud is to provide 
a confining ram pressure. We will see in the next sec- 
tion that accretion also increases the turbulent velocity 
dispersion, implying that the kinetic pressure term also 
increases when accretion is included. The cloud radius is 
determined by a balance between kinetic pressure and a 
combination of gravity, accretion ram pressure, and wind 
recoil pressure. Thermal pressure support is negligible. 

We also note that although the ambient pressure term 
is of the same form as in Paper I, we assume an ob- 
servationally motivated val ue for the am bient pressure, 
^ambA'B = 3x 10 4 K cm" 3 (|McKedll999T ). This includes 
thermal and turbulent pressure but neglects magnetic 
and cosmic ray pressure, since magnetic fields and cos- 
mic rays permeate both the cloud and the ambient ISM. 
We have also adjusted the ambient pressure upwards by 
a factor of two because GMCs form in overdense regions 
of the ISM where the hydrostatic pressure is higher than 
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average. In Paper I, P am b was chosen to be artificially 
high to ensure that the cloud would start its evolution 
in hydrostatic equilibrium. This choice was made to ac- 
count for the weight of the gaseous reservoir that was 
not explicitly included. In practice, by choosing a lower 
value for -P am b, we find that the ambient pressure term 
is subdominant for most of the evolution of the cloud. 
This is expected, since we now correctly account for the 
pressure of the reservoir through the term proportional 
to 7/a- Once accretion halts, the cloud is left out of pres- 
sure equilibrium and must expand to match the ambient 
pressure. This effect is seen most clearly in the second 
column of Figure [3] 

2.3. Energy Equation 

In Appendix [Cj we derive the time evolution equation 
for the total energy of the cloud, 



d£ c \ _ M cl 
dt ~ M cl [tcl 



(1 - r,%)W) 



47rP ambj Rci^ci 



GM cl M cl M cl R cl 
-X 1 



3-kf 



M C \R C \ 



Rci{M ci v cj - £M. 



acc ^csc 



(17) 



3 • 3 • 

+ ^ M accO- r 2 os + -M acc a c} + 7M acc <u esc ) 

— a/M acc i?^ — 3M acc a cl + Q c i - £ cl , 

where Q c \ and C c \ are the rates of energy gain and loss 
due to H II regions and turbulent dissipation, respec- 
tively, £ c \ is the tot al en ergy due to the presence of the 
cloud (see equation |C5j). cr rcs is the velocity dispersion 
of mater ial that is being accrete d, x 1 S given by Equa- 
tion IC81 7 is given by equation (|C15[) , and tp is a free 
parameter that sets the amount of energy available to 
drive accretion driven turbulence. The parameter tp is 
the only adjustable constant in our model that is not 
constrained by the results of simulations or observations 
and must be tuned to reproduce the observed properties 
of clouds. The evolution of the cloud is very sensitive to 
tp and we justify our fiducial choice, tp = 0.75, in §[3] 

This equation governs the global energy budget of the 
cloud. Each term has a straightforward physical expla- 
nation. The term (M c i/M c \)£ c \ is the rate of change of 
the cloud's energy as mass is advectively added to or car- 
ried away from the cloud. Similarly, (M c i/M c i)(l-r?|)W 
is the rate of change of the gravitational and magnetic 
energy due to changes in the mass of the cloud. The next 
term is the rate at which external pressure does compres- 
sional work on the cloud. This is followed by a term that 
accounts for the gravitational work done on the cloud by 
the reservoir as the cloud expands and contracts. The 
following term represents the rate at which mass inflows, 
outflows, and external thermal and turbulent pressure, 
respectively, do compressional work on the cloud. The 
next term, which is proportional to tp, represents the rate 
of kinetic energy injection via stirring of turbulence by 
accreted material. This is followed by two terms that are 
proportional to M acc , which account for the fact that in 
the frame comoving with the motions of material in the 
cloud, accreted material is moving at the transformed 



velocity, v rcs — v, different than the velocity of the reser- 
voir material in the rest frame, v res . Lastly, Q c \ and C c \ 
are the rate of energy injection by star formation and the 
rate at which energy is radiated away, respectively. 

Noting that turbulent motions carry no net radial flux 
of matter and recalling tha t we had set £>turb = 0.67tmb, 
we may evaluate equation (|C5[) and obtain for the total 
cloud energy, 



£o\ = 7 :a I M cl R 2 cl + 2AM cl a 2 cl 



-Mac 1 



§o'(l - V%) 



X 



GMl 



(18) 



Taking the time derivative of this expression, substitut- 
ing into equation (|17[) , and nondimensionalizing as in 
§ 12.21 yields a time evolution equation for a, 



R'R 1 ' 



VG 



Ma 



VA 



MR' 

(Mi?)!/2 C 



VP- 



R 2 R' 
Ma 



M L CC R' 2 



Ma 



(3-1.5tp)M^ c a , <^M ac 



aiM 

+ g_c 

ajMa 
where ? = a les /a xes>0 , 

g-c 



VD- 



Ma 



vi- 



fRa 



(19) 



Rcifi{Qcl — £cl) 



M cl,0CTcl, 



and we define the constants, 



Vd = 



nJ2 
oa rcs,0 

2 °I (7 cl,0 ' 



10 



vi =- 



aia V ir,o 



(20) 

(21) 
(22) 



Here, jjd is propotional to the ratio of the initial turbu- 
lent kinetic energy in the reservoir and the initial turbu- 
lent kinetic energy in the cloud and rjj is proportional to 
the ratio of the initial kinetic energy due to the infall of 
the reservoir to the initial turbulent kinetic energy of the 
cloud. 

Since motions in GMCs are highly supersonic, the 
internal structure of a typical cloud is characterized 
by strong shocks. Because clouds have short cooling 
timescalcs, the shocks present throughout GMCs must be 
radiative. The braking of turbulent motions via radiative 
shocks has be en extensively studied in numerical simu la- 
tions (see e.g. iMac Low etaT]|1998t iStone et al.1ll998D in 
which the turbulent dissipation timescale is found to be 
tdis = Eturb/E = k\i n /a c \ where k is a constant of order 
unity and A;„ is the characteristic length scale of tur- 
bulent energy injection. The simulations of IStone et al.l 
(|1998f ) give k = 0.48 and E tm b = 2AM cl a 2 cV Motivated 
by this r esult and using a scalin g argument given by 
IMatznerl (|200l and IMcKeel (fl989h . we assume that the 
dimcnsionless rate of energy loss is given by 



C 



Ma 3 



R 



(23) 
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Here r\ v is a constant of order unity that depends on the 
nature of MHD turbulence in the cloud and we assume 
0in = A; n /4i? c i. The factor of 4 in our expression for </>i n 
comes from the fact that the largest wavelength mode 
supported by the cloud is A max = 4i? c i, corresponding to 
net expansion or co mpression of the clo ud. We make use 
of the simulations of iStone et al.l (|1998t ) to calibrate this 
expression. For the run most resembling real molecular 
clouds, we find rj v = 1.2. We follow Paper I in adopting 
4>- m = 1.0 below. This is moti vated by the results of 
Brunt et al.1 (12009T) fbut se e also lOssenkopf fe Mac Low! 
2002t iHever fc Brunt! 120041 ) who compared the velocity 
structure of observed clouds, where Ai n cannot be directly 
observed, with the velocity structure of simulated clouds, 
where A; n is known a priori, and found Ai n > R c \. 

Comparing our velocity dispersion evolution equation 
(equation |19| ) to the corresponding equation given in 
Paper I, we see there are four new terms proportional to 
M' &cc . In practice, we find that the primary effect of ac- 
cretion on the energy balance of the cloud is to increase 
the turbulent velocity dispe rsion via the terms propor- 
tional to ip. We will show in § !4.2l that the velocity disper- 
sion is set by a balance between the decay of turbulence 
and energy injected by accretion and star formation. 

2.4. Star Formation and H II Regions 

Star formation is able to influence the evolution of the 
cloud by ejecting mass and by injecting turbulent kinetic 
energy as expanding H II regions merge with and drive 
turbulent motions in the cloud. The first mechanism 
is accounted for in our models by including an ionized 
wind that decreases the mass of the cloud and confines 
the cloud by supplying recoil pressure. The second mech- 
anism is accounted for by the Q c \ term in equation (JTTJ) 
that represents energy injection by H II regions. 

Since we only know the global properties of the 
cloud, we calculate the rate of star formation by mak- 
in g use of a power-law fit to the star formation law 
of iKrumholz fc McKei (|2005[) . Stars form at a low 
efficiency per free-fall time, consistent with observa- 
tions of star formation in nearby molecular clouds 
(jKrumholz fc Tan] l2007f ). Individual star formation 
events occur once a sufficient amount of mass has ac- 
cumulated to form a star cluster. Star cluster masses 
are found by drawing from a cluster mass function ap- 
propriate for a single cloud (see equation [44] of Paper 
I). We then populate t he cluste r with individual stars by 
picking masses from a iKroupal (|2002l ) IMF. If the total 
ionizing luminosity of the newborn star cluster is suffi- 
cient to drive the expansion of an H II region, we begin 
to track the resulting expansion. 

Once a massive star cluster forms, it photoionizes gas 
in its surroundings and drives the expansion of an H II 
region. Paper I tracked the expansion of individual H n 
regions by assuming the analytic self-si milar solu t ion fo r 
H II region expansion worked out by iMatznerl (|2002t ). 
This solution uses the fact that once an H n region has 
expanded beyond the Stromgren radius, most of the mass 
in the H n region volume is in a thin shell of atomic gas 
at a radius r s h from the center of the H II region. The 
ionized gas in the interior of the shell exerts a pressure 
on the surface of the shall, causing the shell to accelerate 
outwards. The shell evolution equation derived from this 
analysis admits a self-similar solution for the expansion 



of the H II region. This self-similar result does a good job 
of predicting the expansion if there is no characteristic 
scale in the problem. 

However, the introduction of radiation pressure leads 
to a characteristic radius, r c h, and time, i c h, at which the 
gas pressure and radiation pressure at the inner surface 
of the shell are equal. Radiation pressure is the dominant 
force driving the expansion of the ionized bubble when 
r s h < r r h and gas pressure dominates when r s h > r c h. 
IKM09I modified the theory of lMatznerl (|2002l ) to account 
for the effect of radiation pressure in the initial stages of 
the expansion. They derived an explicit functional form 
for r c h and i c h in terms of the bolometric and ionizing 
luminosity of the central star cluster, properties of the 
molecular cloud, a nd fund amental constants (see equa- 
tions [4] and [9] in IKM09T ). The numerical value of r c h 
and t c h depends on the the bolometric luminosity of the 
central star cluster and the ionizing photon flux of the 
central star cluster. The value we choose for /trap, a 
factor that accounts for the trapping and reradiation of 
photons as well as the trapping of main sequence winds 
within the neutral shell, and </>, a factor that accounts 
for the absorptio n of ra diation by dust, are the fiducial 
values quoted bv IKM091 

Defining the dimcnsionlcss variables x^ — r s h/r c h and 
7sh = t/t c h , the equation of motion for the shell reduces 
to (|KM09D 



This assumes that gas in the neighborhood of the H II 
region follows a density profile proportional to r _1 , effec- 
tively placing the H II region in the center of the cloud. 
This accounts for the fact that H n regions form in over- 
dense regions of the cloud. In practice, we solve equa- 
tion (|24j) numerically to obtain x s h and dx s h/d,T s h and 
thus f s h(t) and r s h(t). 

In a gas pressure driven H n region, the force exerted 
on the expanding bubble is twice the recoil force pho- 
toionized mate rial imparts on the cloud as it is ejected 
( Matzner 2002) . The photoevaporation rate can then be 

straightforwardly calculated via M c j = — p s \ 1 /2c u . Here 
p s h is the momentum of the shell and p s h is the force 
acting on the shell. In a radiation pressure driven H II 
region, the total force is given by the sum of the gas pres- 
sure and radiation pressure forces, p s h = p gas + Pr&d- At 
early times, when r s h *C f c h, the radiation force domi- 
nates, so p ra( j Pgas- Thus, If we calculate the mass 
ejection rate from the the total force acting on the shell, 
we will overestimate the mass ejection rate in a radiation 
pressure dominated H n region. To correct for this effect, 
we modify the analysis of Paper I by only including the 
gas pressure force when we calculate the mass ejection 
rate, 



M cj = -ff (25) 
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where, 

P gas = (^) 1/2 2.2k B T n rie (26) 

= 3.3 x 10 28 L 39 x s Y 2 dynes. (27) 

Since we do not include the dynamical ejection of ma- 
terial as H II regions break out of the cloud surface, 
this is formally a lower limit on the true mass ejection 
rate. Since cloud s are clumpy and somewhat porous 
(|Lopez et al.l 120 111 ) , we expect to make little error by 
neglecting dynamical ejection. 

Once the stars providing 50% of the total ionizing lu- 
minosity of the central star cluster have left the main 
sequence, the H n region enters an undriven momentum- 
conserving snowplow phase. When the expansion veloc- 
ity of the H II region is comparable to the turbulent ve- 
locity dispersion, we assume that the H II region breaks 
up and contributes turbulent kinetic energy to the cloud. 
H II regions can merge with the cloud during cither the 
driven or undriven phases. If the radius of the H II region 
is greater than the radius of the cloud and the expansion 
velocity of the shell is greater than the cloud escape ve- 
locity, we say the H n region disrupts the cloud and end 
the global evolution. 

If an H II region merges with the cloud at time t = t rn 
when its radius is r s h = r m , the rate of energy injection 
from a single H II region is given by 

/ \ 1/2 

G ci = 1.67] E Ti(t m ) (^-j 5(t - t m ). (28) 

Here 7i(t m ) — p s hC c i/2, t\e paramaterizes the efficiency 
of energy injection, S(t) is the Dirac delta function, 
and the factor of 1.6 arises because magnetic turbu- 
lence is slightly sub-equipartition compared to kinetic 
turbulence. The factor {r m / R c \) x / 2 accounts for the 
more rapid decay of turbulence when the driving scale 
is smaller (see Paper I). 

2.5. Mass Accretion 

Consistent with the analysis in Appendix [B] we treat 
the reservoir as a gravitationally unstable spherical cloud 
undergoing collapse. The cloud is primarily composed of 
atomic gas in both warm and cold phases. We expect 
that as the reservoir collapses, material that is accreted 
onto the cloud will cool and become molecular. We ap- 
proximate that the reservoir has approximately constant 
surface density, S rcs . To find an upper limit on the mass 
accretion rate, we can assume that the reservoir is under- 
going collapse in the limit of zero pressure. It is straight- 
forward to show that the resulting accretion rate onto 
the central condensation is given by 

Maccff = — G 2 SLi 3 (29) 

7T 

Since the estimate for the accretion rate given in equa- 
tion [25] does not take into account pressure support, 
it is lik ely an overestimate of the true accretion rate. 
Indeed, iMcKee fc Tanl (120031 ) considered the inside out 
collapse of equilibrium polytropic molecular cloud cores 
and, in the case of k p = 1, found the same scaling with 



time and surface den s ity bu t a substantially lower coeffi- 
cient. iTan fc McKeej (|2004D argued that subsonic inflow 
was a more realistic initial condit ion condition than th e 
static equilibr ium solution used bv lMcKee fc Tanl (|2003|) . 
IHunterl (|1977f ) found that the set of solutions for the col- 
lapse of an isothermal sphere starting at an infinite time 
in the past. The solution with an infall velocity of about 
c c i/3 corresponds well to the results of the simul a tion o f 
the formation o f a pr imordial star by I Abel et al.l (|2002f) . 
ITan fc McKeel (|2004l ) adopted this solution, noting that 
it has an accretion rate that is 2.6 time s greater than that 
for a static initial condition (lShulll977T ) when expressed in 
dimensionlcss form. Our problem is quite different from 
Hunter's, since an equ ilibrium density gra dient k p = 1 
corresponds to 7 = (|McKee fc Tanll2003f ) rather than 
Hunter's 7 = 1. Nonetheless, we assume that the accre- 
tion rate for our problem is also 2.6 times greater than 
that for a static initial condition and find 

Macc,TM04 = 10.9G 2 £L* 3 - (30) 

Clearly, this result is uncertain, and magnetic fields in- 
troduce further uncertainty. Fortunately, we find that 
varying the numerical coefficient in equation [30] does not 
affect the qualitative nature of the results discussed be- 
low. 

To model the effect of a finite gas supply, once the 
total mass of gas that has fallen onto the cloud ex- 
ceeds the total mass of the reservoir, M res , we set 
M acc = 0. When comparing with galactic populations 
of GMCs, we set M Tes = 6 x 10 6 M B , the observed 
upper mass cutoff for GM Cs (| Williams fc McKed II9971: 
iFukui fc Kawamural[2010l ). This may underestimate the 
true upper mass since fragmentation may lead to a range 
of reservoir masses. Since 6 x 10 6 M Q GMCs are rela- 
tively rare, we set M ros = 2 x 10 6 Mq for the runs pre- 
sented in Figure [3] and Table [5] as ~ 10 6 M Q is a more 
typical GMC mass. 

If we also assume that the atomic reservoir is virialized 
such that the virial parameter of the reservoir, a res , is 
constant with radius, we find 

a res = 0AaUsGX Ies t (31) 

Fiducially, we take a ros = 2.0, corresponding to a 
marginally gravitationally bound reservoir. This param- 
eterization assumes that the reservoir satisfies an internal 
linewidth-size relation of the form cr re s(f) oc r 1 / 2 . The 
increase of the reservoir velocity dispersion with time re- 
flects that material originates at increasingly larger radii. 

The precise normalization of the mass accretion law is 
a major source of uncertainty in our modeling. Since the 
length scales over which material is swept up into the 
cloud t hrough the reservoir approach g alactic dynamical 
scales (jDobbs et al.ll201ll : lTaskerll2"011f L a more complete 
treatment would require tracking the gas dynamics from 
the scale of an entire galaxy down to the the scale of the 
reservoir. This would also allow us to self-consistently 
model the end of accretion onto the cloud instead of as- 
suming that mass accretion cuts off abruptly. In a forth- 
coming paper, we plan to add our molecular cloud mod- 
els to a simulation of gas dynamics in a galactic disk to 
model the gas reservoir for a population of GMCs. 

Although the normalization of the accretion law is 
somewhat uncertain, we can make use of observations 
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of the gas content of nearby galaxies to estimate E res , 
the surface density of gas in the reservoir feeding the 
cloud. The average H I surface density in the inner disk 
of the Milky Way is observed to be approxim ately con- 
stant, ~8M Q pc~ 2 (|Kalberla fc Dedesl 120081) . Beyond 
a galactocentric radius of ~ 10 - 15 kpc, the H I surface 
density exponentially decreases with radius, however, few 
GMC s are observed in the outer Milky Way (|Hever et al.l 
I2001D or beyond the o ptical radius of nearby galaxies 
(|Engargiola et al.ll2003t Bigiel et al.l |2008(). Similar satu- 
rated mean H I surface densities are observed in nearby 
galaxies, except in the central regions of some galaxies 
where the ISM becomes ful ly molecular and t he H I sur- 
face density goes to zero (|Lerov et al.l 120081 ) . For this 
reason, we adopt E ros = 8 M pc~ 2 as a fiducial atomic 
reservoir surface density typical of the bulk ISM of local 
star forming galaxies. 

Although the mean atomic surface density in nearby 
galaxies is as low as 8 M Q pc -2 , the atomic ISM is ob- 
served to be clumpy, with overdense regions reaching sig- 
nificantly higher surface densities. These r egions may be 
assoc iated with spiral arms, as in M 33 (jThilker et al.l 
EOOl), or_dr iven by gr avitational instability, as in the 
LMC (lYang et al.H200W) . For this reason, we also explore 
the behavior of molecular clouds accreting from higher 
surface density gas, S ros = 16 M Q pc -2 . Since massive 
molecular clouds are universally obs erved to be associ- 
ated with high s urface density gas (jWong et al.l 120091 : 
llmara et al.l 120 111) , we expect there to be marked differ- 
ences between molecular clouds that accrete from high 
surface density gas and clouds that accrete from low 
surface density gas. Although the gas will be primar- 
ily atomic at a gas surface density of 16 M© pc~ 2 , we 
expect that there should be some diffuse "dark" m olecu- 
lar gas (jKrumholz et al.ll200l IWolfire et al.ll20Tot) . Thus 
the reservoir is not necessarily completely atomic, but 
instead primarily composed of atomic gas. 

2.6. Numerical Scheme 

Equations Q, (JTUJ) , and (fT9")) constitute a system of 
coupled, stochastic, nonlinear ordinary differential equa- 
tions in M, <j, and R. We solve these equations by using 
a straightforward Euler integration with an adaptive step 
size. The precise order in which we update cloud proper- 
ties is as follows. After calculating the instantaneous star 
formation rate, we calculate M' by summing the compo- 
nents due to ionized winds and mass accretion. Next, 
we calculate the rate of turbulent dissipation using equa- 
tion (j2"B"]) . We then calculate the ratio of the mas s 
doubling time to the free-fall time, using equation (|B5|) . 
We use £ to calculate a', /, and £ by interpolating on 
precomputed tables. Since a' depends on R" , we first 
evaluate R" using equation (|TU)) and then compute a' 
using equation (fT9"]) . Next, we check if R, M, or a will 
change by more that 0.1 % using the current value of the 
time step. If we detect a change larger than this, the time 
step is iteratively recalculated using a new time step half 
the size of the original until the fractional changes in i?, 
M, and a are smaller than 0.1%. Next, we calculate 
R, M, and a at the new time step, update the state of 
any H n regions created in previous time steps, and then 
create new H II regions using the procedure described in 
§ 12.41 If the time step did not need to be reduced, we 
increase the size of the time step by 10%. 



Cloud evolution can be terminated if one of three con- 
ditions is satisfied: 

• The time step is less than 10 -8 of the current evo- 
lution time (i.e. Ar/r < 10~ 8 ). 

• The mean visual extinction falls below Av,min — 
1.4, corre sponding to the CO dissociati on threshold 
found bv Ivan Dishoeck fc Blackl (| 19881) . 

• An H II region envelops and unbinds the cloud, i.e. 
if r sh > R c i and f sh > 

We use the phrases collapse, dissociation, and disruption, 
respectively to describe these scenarios. The dissociation 
thres hold depends on th e ambient radiation field. How- 
ever, [W<5fire]iLlI] (120101 ) found that the CO dissociation 
threshold varies by only a factor of two when the intensity 
of the ambient radiation field varies by an order of mag- 
nitude, so we neglect variations in the radiation field and 
for consistency with Paper I, adopt Ay, m in = 1.4. The 
surface density corresponding to the dissociation thresh- 
old depends on the assumed dust to gas ratio and thus 
on the metallicity. For solar metallicity, a surface density 
of 1 g cm~ 2 corresponds to Ay = 214.3. Since we use a 
dissociation threshold based on CO rather than H2 to 
define the end of the cloud's life, we may miss the fur- 
ther evolution of a diffuse molecular cloud where most of 
the carbon is neutral or singly ionized but the hydrogen 
is still molecular. 

These halting conditions probably oversimplify the 
true end of a cloud's evolution due to our assumption 
of spherical symmetry and homology. For the case of 
collapse, it is more likely that the cloud would undergo 
runaway fragmentation rather than monolithic collapse. 
In the case of dissociation, even if the mean surface den- 
sity drops below the point where CO can no longer re- 
main molecular, that does not preclude the possibility 
that overdense clumps might retain significant amounts 
of CO. Finally, for the case of disruption, even if an H II 
region delivered a large enough impulse to unbind the 
cloud, it may simply be displaced as a whole, or be dis- 
rupted into multiple pieces which would then evolve in- 
dependently. It is also likely that if the cloud is disrupted 
while still actively accreting, the cloud would inevitably 
recollapse since it is unlikely that an H II region would 
have enough kinetic energy to unbind the reservoir. 

Since we must necessarily use simple criteria to halt 
the cloud evolution, our estimates of cloud lifetimes pre- 
sented below arc strictly lower limits to the true lifetime 
of a cloud. Both the disruption and dissociation criteria 
do not preclude the presence of overdense clumps that 
may survive the destruction events. However, our esti- 
mates of cloud lifetimes are appropriate for the lifetime 
of a single monolithic cloud. Any overdense clumps that 
do survive would represent entirely different clouds that 
would evolve independently of each other. 

2.7. Input Parameters 

To complete our model, we must choose a set of param- 
eters and initial conditions to fully determine our cloud 
evolution equations. In Table [T] we have listed the vari- 
ous fiducial parameters we have chosen for our model as 
well as the references from which we derive our choices. 
Some of these parameters are motivated by observations, 
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Table 1 

Fiducial Parameters 



P HI' 3,1110 1 ( 'I' 


Value 


.H cfcrcncc 


«vir,0 


2.0 


Blitz et al. (2007) 


^cl,0 


60 M pc~ 2 


Hever et al. (2009) 


c s a 


0.19 km s -1 




c„ b 


9.74 km s" 1 


McKee & Williams (19971 


IB 


0.5 


Krumholz & McKee 755o5~) 


IE 


1.0 


Paper I 


Tv 


1.2 


Paper I 


<t>in 


1.0 


Brunt et al. (2009) 


f 


0.75 


This work 


A 

,min 


1.4 M Q pc~ 2 


van Dishoeck & Black (1988) 


Ores 


2.0 






3 x 10 4 K cm" 3 


McKee (1999) 


M ch0 


5 x 10 4 M Q 





a Assumes T = 10 K and /x = 2.3 
b Assumes T = 7000 K and fi = 0.61 

others by the results of simulations, and one parameter 
(tp, sec §[3]) is left free. 

The cloud initial conditions can be computed given 
an initial cloud mass along with our ass ume d value for 
own and £ cl , . iMatzner fc McKed (|2000D and lFall et all 
(|2010t ) found that protostellar outflows are energetically 
important when M c \ < 10 4 5 M Q . Above this mass, they 
contribute negligibly. Thus, if we choose a low initial 
mass, we may be underestimating the amount of turbu- 
lent energy injection by star formation feedback at early 
times since we do not account for protostellar outflows. 
For this reason, we choose a relatively large initial mass, 
M c \,o = 5 x 10 4 M Q . For reference, given our choice 
of initial virial parameter and surface density, this cor- 
responds to, i? c i,o = 11.5 pc, a c \fl = 2.7 km s" 1 . and 
t C x,o = 4.1 Myr. While this is larger than some local 
molecular clouds like Taurus or Perseus, it is still much 
smaller than the mass of the molecular clouds where 
most of the star-formation in local galaxies occurs. This 
choice also ensures that we are accounting for the bulk 
of the energy available from star formation feedback due 
to photoionization, ionized gas pressure, and radiation 
pressure. 

3. MODELS WITH ACCRETION ONLY 

Before beginning full simulations using our method, we 
must first test the behavior of our model as we vary the 
free parameter tp introduced in the derivation of the en- 
ergy evolution equation. For this purpose, we have run 
our model with star formation feedback disabled. The 
only physical mechanisms modeled in these tests are ac- 
cretion and the decay of turbulence. Since there is no 
random drawing from the stellar or cluster IMF in these 
runs, the results are fully deterministic. These simplified 
models allow us to understand how the results depend on 
our choice for the tunable parameter tp and provide phys- 
ical insight that will be useful in interpreting the results 
of the more complex and stochastic runs that include 
feedback. 

The energetics and virial balance of our cloud models 
depend critically on the parameter tp. Broadly speak- 
ing, ip controls the amount of turbulent kinetic energy 
injected by the accretion flow. For the case <p = the ac- 
cretion flow contributes the minimum possible amount of 
turbulent kinetic energy. This means that accreted mate- 
rial cannot contribute significantly to turbulent pressure 
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Figure 2. Cloud surface densities (bottom row), virial parame- 
ters (second row), velocity dispersions (third row), and radii (top 
row) for 400 different runs, each with a different choices for tp, as 
indicated in the color bar. Star formation was turned off for all 
runs. 

support, since the accreted material is maximally sub- 
virial. With this choice, as the cloud accretes mass, its 
energy budget must become increasingly dominated by 
self-gravity. Once this happens, internal pressure sup- 
port is negligible and the cloud must inevitably undergo 
gravitational collapse. 

Alternatively, we could set ip > 0. If tp is small, the 
turbulent kinetic energy of a newly accreted parcel of 
gas would still be small compared to the gravitational 
potential energy of the gas parcel. Thus, once the cloud 
is primarily composed accreted gas, the cloud will un- 
dergo gravitational collapse, although on a slightly longer 
timescalc than in the tp = case. At some larger value of 
(f, accretion contributes a net positive amount of energy 
to the cloud, balancing out the negative gravitational po- 
tential energy of the newly accreted material. The could 
will still collapse with this choice since turbulent motions 
quickly decay away. As tp is increased further, we should 
eventually find that at some critical value, <p = tp cr it, ac- 
cretion drives turbulent motions with sufficient vigor to 
avoid the gravitational collapse of the cloud entirely. 

The results of test runs with different choices of tp are 
presented in Figure [5] The time evolution of the cloud 
surface density, virial parameter, velocity dispersion, and 
cloud radius is plotted for a selection of clouds evolved 
with different choices for tp. Each line depicts the time 
evolution of a cloud property and is color coded by the 
value of tp chosen. It is obvious that the value of tp can 
strongly influence the resulting evolution. 

If tp = 0, the cloud experiences global collapse in a 
free-fall time. Initially, the cloud velocity dispersion 
decreases, but inevitably the gravitational term in the 
velocity dispersion evolution equation, proportional to 
— R' / R 2 , becomes dominant, and the velocity dispersion 
begins to diverge. The fact that a c \ diverges as R c \ goes 
to zero is an artifact. In reality, the highest density re- 
gions would independently fragment and collapse and the 
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cloud would never undergo a monolithic collapse. 

As we increase ip, the cloud is able to support itself 
against collapse for longer periods. Near <p = 0.5, accre- 
tion brings in net positive energy but turbulent dissipa- 
tion wins out, and the cloud still eventually collapses. At 
a critical value, p> cr it — 0.8, accretion driven turbulence 
alone is sufficient to hold up the cloud against collapse for 
as long as the reservoir continues to supply mass to the 
cloud. The mass, radius, and velocity dispersion of the 
cloud increase in such a way as to maintain a constant 
virial parameter and surface density. 

Since we expect that gas motions driven by accret- 
ing dense clumps should be at least somewhat correlated 
with the motions of the infalling clumps, we do not ex- 
pect a physically realistic choice of ip to be very close to 
zero. On the other hand, a model in which a cloud is en- 
tirely supported by accretion driven turbulence seems to 
preclude the possibility that a significant fraction of the 
kinetic energy of infalling gas is radiated away in an ac- 
cretion shock. For this reason, we rule out as unphysical 
runs with ip sa and ip > <p cr it ■ The precise value of ip we 
will use in our models that include star formation below 
depends on uncertain details of the accretion and mix- 
ing of infalling gas. In practice, we find that even with 
the energy provided by star formation feedback, clouds 
generally undergo free-fall collapse or reach unreasonably 
high mean surface densities once they are primarily com- 
posed of accreted material if we choose <p < 0.7. Since 
clouds are generally not observed to be in global free- 
fall collapse, we instead pick a value somewhat higher 
that this, p = 0.75 for our fiducial models. This splits 
the difference between accretion contributing a negligi- 
ble amount of energy to the cloud when ip = 0.5 and 
accretion contributing the maximum possible amount of 
energy when tp = 1. We will see below that our fidu- 
cial choice broadly reproduces the observed properties of 
molecular clouds in the Milky Way and nearby galaxies. 

4. MODELS WITH ACCRETION AND STAR FORMATION 

Feedback by the action of ionizing radiation emitted 
by newborn stellar associations alters the evolution of a 
GMC after the birth of the first massive star cluster. The 
source of energy provided by massive star formation can 
be a significant component of the energy budget of the 
entire cloud. For the remainder of this paper, we consider 
models with the star formation prescription described in 
§EH turned on. 

4.1. Overview of Results 

We have run two sets simulations with parameters cho- 
sen to model conditions in interarm (E res = 8 M Q pc~ 2 ) 
and spiral arm (£ rcs = 16 M pc~ 2 ) regions. Besides 
the two different choices for the ambient surface density, 
all other parameters and initial conditions are identical. 
The time evolution of a subsample of runs are plotted 
in Figure [3] and average properties of the full sample are 
presented in Table [5J 

The most striking result of our comparison is that the 
final mass of our model molecular clouds depends on the 
assumed mass accretion history. Clouds evolved with a 
low accretion rate, corresponding to conditions in inter- 
arm regions, grow larger than 10 5 M Q less than 30% of 
the time and very rarely reach masses comparable to the 



most massive GMCs in the local group. The vast major- 
ity of clouds are instead disrupted by an energetic H n 
region within a few crossing times. The clouds attain 
a quasi-equilibrium configuration in which mass accre- 
tion is roughly balanced by mass ejection. Clouds avoid 
global collapse by extracting energy from the expansion 
of H II regions. 

The evolution of the clouds is characterized by discrete 
energy injection events due to the formation of a sin- 
gle massive star cluster. Once a cluster forms, it ejects 
a wind and launches an H n region. The recoil force 
of launching the wind leads to an overall confining ram 
pressure, causing the radius to decrease and the surface 
density to increase. Once the star cluster burns out, the 
H II region expansion decelerates and then stalls. When 
the expansion velocity of the H n region is comparable 
to the cloud velocity dispersion, the kinetic energy of the 
expanding H n region is converted into turbulent kinetic 
energy, causing a spike in the turbulent velocity disper- 
sion. The turbulent kinetic energy exponentially decays 
away over a crossing time, but the temporarily elevated 
velocity dispersion increases the turbulent kinetic pres- 
sure, causing the cloud to expand. This leads to oscilla- 
tions in the cloud radius and mean surface density. On 
the whole, clouds that are not quickly disrupted by H II 
regions, are able to survive as quasi- virialized objects for 
several crossing times before they are cither disrupted or 
dissociated. 

Clouds evolved with a higher ambient surface density, 
typical of spiral arm regions in the Milky Way, exhibit 
significantly different behavior. Since these clouds ac- 
crete mass much faster than in the low surface density 
runs, they are not able to attain steady state between ac- 
cretion and ejection of mass. While some clouds are still 
destroyed by energetic H II regions early in their evolu- 
tion, over 90% of these clouds were able to accrete their 
entire reservoir after 25 Myr. At this point, the clouds 
are generally quite massive, ~ 1.5 x 10 6 Mq. Once ac- 
cretion is shut off, the clouds are no longer confined by 
accretion ram pressure and lose a portion of the power 
that had been driving turbulence. For this reason, the 
velocity dispersion decreases in response to the loss of ac- 
cretion driven turbulence, and the cloud radius expands 
in response to the loss of the confining pressure provided 
by accretion. Before the cloud can dissociate, it attains 
pressure balance with the ambient ISM at a lower ve- 
locity dispersion and larger radius. For the next 20 to 
30 Myr, the clouds evolve in much the same way as the 
massive cloud models considered in Paper I. The clouds 
can be supported against self-gravity for many dynam- 
ical times by forming stars and launching H n regions. 
Particularly energetic H n regions can disrupt the clouds 
and excursions to low surface density can dissociate the 
clouds. The lifetime of these clouds is thus set by the 
amount of time they can accrete. This may imply that 
spiral arm passage times set GMC lifetimes, although 
further work is needed to clarify this tentative conclu- 
sion. 

The star forming properties of the two sets of models 
are also somewhat different. In the interarm case, the 
star formation efficiency, 
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Average Properties of Model Clouds 



^res (^life) 

[M Q pc- 2 ] [Myr] 


<M max ) 
[Mq] 


(-^ejected) 

[Mq] 


(e) (6ff) 
[%] [%] 


-^dissoc ^disrupt 


8 26.2 ±29.8 
16 52.6 ±16.8 


(3.7 ±4.9) x 10 s 
(1.3 ±0.4) x 10 6 


(3.4 ± 5.3) x 10 5 
(1.2 ±0.4) x 10 6 
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Figure 3. Cloud surface densities (bottom row), star formation rates (second row), virial parameters (third row), velocity dispersions 
(fourth row), and radii (top row) for a set of 40 clouds. Half of the clouds were evolved with a low surface density characteristic of the 
bulk of the atomic ISM and the other half were evolved with a high surface density, characteristic of overdense regions in the ISM. The 
two different choices of E res are marked at top. Accretion was shut off once 10 6 Mq of material had been processed through the accretion 
flow. The ambient surface density, and thus the accretion history, strongly affects the resulting cloud evolution. 
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is only 5% while in the high surface density case, e is 
somewhat larger, approximately 8%. This can be entirely 
attributed to the difference in lifetimes between the two 
sets of models. In the low surface density case, most 
clouds are only able to survive one or two crossing times 
and thus can only convert a small fraction of their mass 
into stars before they are dissociated or disrupted. The 
clouds evolved with a high ambient surface density are 
able to survive for many crossing times and convert a 
larger fraction of their gas into stars. An even larger 
fraction is ejected via photoionization. However, for both 
models, the star formation efficiency per free-fall time, 



M* 
Mats 



(33) 



is low, around 2%. This is not surprising, as a low star 
formation efficiency per free-fall time is one of the basic 
assumptions of our model. 

4.2. Energetics of Star Formation Feedback Versus 
Mass Accretion 

GMCs exhibit highly supersonic turbulence. There is 
no agreement in the literature about what drives these 
motions, which numerical models of compressible MHD 
turbulence indicate should decay if left undriven. Some 
authors suggest that the primary energy injection mech- 
anism is some sort of internal star for mation feedback 
process, such as proto stellar outflows dLi fe Nakamura 
2006 i lWang et al.l2010fl. expanding H H regions (iMatznei 

3 k 



2002 ) . or supernovae fMac Low fc Messed 120041 ) . Oth 
ers sugg est that turbulence is driven externally via mass 
inflows (jKlessen fc Hennebelld [2010). Comparing the 
amount of energ y injected by diffe rent forms of star for- 
mation feedback, iFall et all (|2010f ) found that at typical 
GMC column densities, the dominant stellar feedback 
mechanism is H II regions driven by the intense radia- 
tion fields emitted by massive star clusters. Using our 
models, we can compare the importance of accretion rel- 
ative to H II regions in the energy budget of GMCs. 

To find the total energy injected by accretion, we make 
use of our knowledge of the total energy of the cloud as 
a function of time. At the end of time step j we use 
equation (|18p to calculate both the total cloud energy, 
S c \j , as well as what the cloud energy would have been if 



we had set M a 
difference, 



for that time step, £, 



cl|M aoc =0- 



cl,j 



lM acc =0 



The 



(34) 



is the total energy added by accretion during that time 
step. The total energy injected by accretion over the 
cloud's lifetime is just the sum of the contributions of 
each time step, 

£&cc — ^ ^ £&cc,j ■ (^^) 

i 

The energy injected by H II region i, fmi.i, can be 
found by integrating the rate of energy injection by a 
single H II region with respect to time. This is, 



1/2 



(36) 



where r m j is the radius of H II region i when it merges 
with the parent cloud and i? c i , is the radius of the cloud 



as a whole when H II region i merged with the cloud. 
To find the total energy injected by H II regions over the 
cloud's lifetime, we simply sum up the contributions due 
to individual H II regions, 



(37) 



The ratio |£h n/£acc| indicates the relative importance 
of star formation feedback to accretion driven turbulence 
to the global energy budget of the cloud. If |£h n/£acc| < 
1, accretion dominates the energy injection; similarly if 
\£n n/£acc| > 1, star formation feedback is the primary 
driver of turbulence. 

The results of this comparison are plotted for both 
choices of the ambient surface density in Figure 2] We 
find that H n regions and accretion contribute approxi- 
mately equal amounts of energy in the low surface den- 
sity runs, while accretion dominates in the high surface 
density runs. In the low surface density runs, stochastic 
effects can be important, particularly for clouds that do 
not last much longer than a crossing time. Thus, in some 
runs, star formation feedback can contribute significantly 
more energy than accretion, while in others star forma- 
tion feedback is negligible. In the runs evolved with a 
high ambient surface density, star formation feedback is 
subdominant, although not completely negligible, in the 
vast majority of runs. 

It is worth pointing out that this result depends on 
the precise value of ip we choose to evolve the clouds 
with. If ip is lower, accretion contributes less energy, and 
star formation can dominate the energy budget. If ip 
is higher, star formation becomes completely negligible, 
and the amount of kinetic energy injection is controlled 
by the mass accretion rate. Since clouds collapse when we 
choose ip much lower than our fiducial value, and shocks 
in molecular gas tend to be strongly dissipative, we do 
not expect the 'true' value of ip to be much different than 
our fiducial value. We thus conclude that one of three 
cases must hold. Star formation may be dominant, but 
only marginally so. Accretion may also be dominant, but 
again, only marginally. It is also possible that star for- 
mation and accretion contribute roughly equal amounts 
of energy. In all three cases, neither star formation or 
accretion is truly negligible. 

5. OBSERVATIONAL COMPARISONS 

5.1. Larson's Laws 

Giant molecular clouds are observed t o obey three 
scaling relations, known as Larson ' s law s (|Larsonlll981l ; 
iSolomon et"aT1ll987t iBolatto et aHl2008[ ). In their sim- 
plest form, Larson's laws state: 

• The velocity dispersion scales with a power of the 
size of the cloud. Subsequent observations have 
shown that this power is about 0.5 (o~ c \ oc i?J?j 5 ). 

• The mass of the cloud scales with the square of the 
radius (constant £ c i). 

• Clouds are in approximate virial equilibrium (a v i r 
of order unity). 

These laws are not independent; any two imply the 
other. At a minimum, an acceptable theoretical model 
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Figure 4. The number of clouds plotted as a function of n/£acc|- In regions of low ambient surface density, accretion and star 
formation are in equipartition, while in regions of high ambient surface density, accretion dominates the energy budget. 



for GMCs should agree with both the scaling and the 
normalization of the Larson scaling relations observed in 
real clouds. We have already seen that clouds maintain 
approximate virial equilibrium as well as roughly con- 
stant surface densities, but we have yet to see whether 
the normalization of the Larson scaling relations for our 
models agrees with the observed Larson scaling relations. 

5.1.1. Equilibrium Surface Densities 

GMCs, both in the Milky Way (jLarsonl 119811: 
ISolomon et all [1987D . and in nearby ex ternal galax- 
ies dBlitz et al 1 l2007t iBolatto et afl l2008h exhibit sur- 
prisingly little varia tion in surface density. For the 
ISolomon et all (|1987l ) sample of Milky Way clouds, this 
was found to be (S c i) = 170 M Q pc -2 . More recent and 
sensitive observations find lower values, closer to (£ c i) = 
50 M Q pc~ 2 in the Milky Way (iHever et all 1 20091) and 
in the Large Magellanic Cloud ([Hughes et al. T2010h . al- 
though these latter estimates depend on a highly uncer- 
tain correction for non-LTE line excitation and the CO 
to H2 conversion factor, respectivel y. Using heterogenous 
data from several nearby galaxies, IBolatto et al.l (|2008l ) 
attempted to extract cloud properties in a uniform man- 
ner and found a typical surface density of 85 M Q pc~ 2 
but with significant variation from galaxy to galaxy. 

Variations in the mean GMC surface density are seen 
when comparing samples from different galaxies. How- 
ever, within a si ngle galaxy there is little variation 
(|Blitz et al.l I2007D . These variations are usually at- 
tributed to differences i n the CO-to-H? conv ersion factor 
from galaxy to galaxy (jBolatto et al.l [2008h . a quantity 
which may d epend on metallicity and the interstellar ra- 
diation field ([Glover et al.ll201Cif) as well as variations in 
turbulent pressure and radiation field in the ambient in- 
terstellar medium. In our runs, we also recover roughly 
constant surface densities (see the second row from the 



bottom of Figured]). 

In Figure El we have reproduced a figure from 
iBlitz et all ([20071 ) that depicts observational results for 
CO luminosities and cloud radii for a sample of clouds 
in the outer Milky Way as well as from several samples 
of extragalactic GMCs. To compare against this com- 
pendium of results, we calculate CO luminosities for our 
model clouds by assuming a constant CO-to-H2 conver- 
sion factor, 



^co 



M cl it I -12 

K km s pc 

3.8 M P 



(38) 



as in IRosolowskv &; Lerovl ([20061 ). This formula ac- 
counts for the presence of helium and assumes a 
constant H2 to CO conversion factor, Xqo — 4 x 
10 20 cm -2 ( K km s -1 ) -1 , twice the value derived for 
molecular clouds within the So lar circle using observa- 
tions of gamma- ray emission ([Strong fc Mattoxl 119961 : 
lAbdo et "all 120101). We cho ose this value to be consis- 
tent with " IBlitz et al.l ([20071 ). who find, using this value 
of Xqo, that all of the GMCs in their sample have virial 
masses comparable to the masses implied by their CO 
luminosity to within a factor of two. 

With our fiducial initial conditions, model clouds in our 
sample begin their lives in the bottom left hand corner 
of Figure at Rci ~ 10 pc. As they accrete and ex- 
pand, clouds move towards the upper right hand corner. 
Clouds end their evolution either through disruption by 
a single H 11 region or by passing below the molecular 
dissociation threshold, indicated by a dashed line in Fig- 
ured Offsets in the distribution of column densities from 
galaxy to galaxy and from the simulated clouds can be 
attributed to variations in Xqo an d uncertai nty in iden- 
tifyin g a unique radius for observed clouds (IBlitz et ah 
2007 ) that have nonzero obliquity (B crtoldi fc McKee 



1992f ). Accounting for variations in Xqo, there is striking 
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Figure 5. Cloud CO luminosity plotted as a function of cloud 
radius. CO luminosities are found by assuming Xqo = 4 x 
10 20 cm -2 ( K km s -1 ) -1 . Solid lines of constant surface density 
arc plotted for 10, 100, and 1000 M Q pc -2 for reference. The 
dashed line of constant surface density corresponds to our assumed 
dissociation threshold. The outputs from a set of 2000 runs were 
used, with £ rea = 8 and 16 M pc -2 and M ICS = 6 X 10 6 M . 
Colors indicate the amount of time model clouds tend to occupy 
a position in Leo — ^cl parameter space. Symbols denote ob- 
served CO luminosities and cloud radi i for galact i c (poi nts) and 
cxtragalactic (open shapes) GMCs. See lBlitz et~aH (120071 ) and ref- 
erences therein for details of the observations. 

agreement between the observed distribution of molecu- 
lar clouds and our sample of simulated clouds. 

The models exhibit a kink in their evolution when 
the reservoir is exhausted and accretion is shut off. 
For this reason, there arc no clouds with Leo > 
10 5 6 K km s _1 pc 2 . Once accretion is shut off, the 
clouds decrease in mass for the remainder of their evo- 
lution. This kink is somewhat artificial since we have 
assumed a fixed reservoir mass and a smooth accretion 
history. A more sophisticated model for the reservoir in- 
cluding a range of reservoir masses would exhibit a con- 
tinuous spectrum of kinks, broadening the region of pa- 
rameter space explored by the models, particularly for 
L CO < 10 4 - 5 K km s" 1 pc 2 . 

The models also exhibit two distinct favored strips of 
parameter space along which they tend to evolve. This 
corresponds to the two different equilibrium column den- 
sities picked out by the two different choices of E res . This 
behavior is clearly seen in the second panel from the bot- 
tom of Figure [H The fact that S c i is sensitive to S res 
follows from dimensional analysis. 

5.1.2. Linewidth-Size Relation 

We next compare our simulated clouds with the 
linewidth-size relation observed t o hold among GMCs as 
a population (jBolatto et al.ll2008h . We are able to repro- 
duce the power law, scatter, and the rough normalization 
in the observed linewidth-size relation. This conclusion 
is unsurprising, since we have already seen that our simu- 
lated clouds maintain roughly constant virial parameters 
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Figure 6. Cloud velocity dispersion plotted as a function of cloud 
radius. The dash-dotted line i s the galactic linewidth-size relation 
found bv lSolomon et al.l l fl987fl . <r v = 0.72R ' 5 . Symbols and color 
coding are the same as in Figure [5] 

and surface densities as they evolve. It is worth noting 
that for our simplified model for the environment of a 
CMC, that the linewidth-size relation corresponds to an 
age sequence. Clouds that live towards the left-hand side 
of the diagram are younger than clouds that live towards 
the right. It is possible that this conclusion is an arti- 
fact of choosing a single reservoir mass. Clouds accreting 
from a population of reservoirs with a continuous spec- 
trum of masses may blur this effect somewhat. We plan 
to revisit this in future work in which we will model the 
global ISM of a galaxy simultaneously with the evolution 
of a population of GMCs. 

There is a small offset when comparing the locus of 
extragalactic and outer Milky Way clouds with our mod- 
els, although there is good a greement betwe e n our mod- 
els and the scaling found by iSolomon et al.l (|1987t ) . For 
a subset of the observational sample, particularly the 
SMC clouds, it is possible that the metallicity of the 
gas in the clouds is so low that CO is no longer a 
good tracer of the bulk of the molecular gas dLerov et al.l 
l200l . Since our models assume perfect sphericity, and 
the observed radius of a prolate or oblate spheroid will 
alway s be smaller than the co rresponding spherical ra- 
dius (jBertoldi fc McKeel 119921 ) . it is also possible that 
the radii predicted by our models overpredict the corre- 
sponding observed cloud radius by 0.1 or 0.2 dex. Lastly, 
it could be the that we overpredict the various pres- 
sures due to photoionization and accretion by assuming 
spherical symmetry. In reality, the wind and accretion 
ram pressure may not necessarily be perfectly spherically 
symmetric, leading to a reduction in the overall confining 
pressure and an increase in the radius. 

5.2. Evolutionary Classification 

The Large Magellanic Cloud is home to one of the best- 
studied samples of GMCs in any galaxy. The LMC's 
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disk-like geometry and face-on orientation offers little 
ambiguity in distance measurements, with t he most accu - 
rate measurements giving ^lmc = 50.1 kpc (|Alvesll2004D . 
A large quantity of high-quality multiwavclcngth data 
has been obtained for the entire disk of the galaxy. In 
particular, the NANTEN 12 CO (J = 1 -> 0) surveys 
and high-resolution f ollowup from the M AGMA 12 CO 
(J = 1 — > 0) survey (jHughes et al.ll20i0t ) have mapped 
the molecular content of the entire disk of the LMC and 
identified 272 clouds that together contain 5 x 10 7 M© 
of molecular gas. When combined with multiwavclcngth 
archival observations of star formation indicators, these 
CO data constitute a snapshot in the evolution and star 
fo rmation history o f a pop ulation of GMCs. 

iKawamura eTaTI (|200l used the NANTEN CO J = 
(1 — > 0) data, along with com plementary Ha photom- 
etry (|Kennicutt fc Hodgelll986|). radio continuum maps 
at 1.4. 4.8. and 8.6 Ghz (toickel et al.|[2005l: iHughes et~aTl 
[2001 and a map of youn g (< 10 Myr) clu sters extracted 
from UBV photometry (jBica et al.l fl996) to investigate 
the ongoing star formation within GMCs in the LMC. 
These authors found a strong tendency for H n regions 
and young clusters to be spatially correlated with GMCs. 
Using this association, the GMCs in their sample were 
separated into three types. Type 1 GMCs are defined to 
be starless in the sense that they are not associated with 
detectable H n regions or young clusters, Type 2 GMCs 
are associated with H n regions, but not young clusters 
in the cluster catalog, and Type 3 GMCs are associated 
with both H II regions and young clusters. 24% of the 
NANTEN sample were classified as Type 1, 50% as Type 
2, and 26% as Type 3. 

Assuming that GMCs and clusters are formed in steady 
state and assuming that young clusters not associated 
with GMCs are associated with GMCs that have dis- 
sipated, one can infer from the NANTEN population 
statistics that GMCs spend 6 Myr in the Type 1 phase, 
13 Myr in the Type 2 phase, 7 Myr in the Type 3 phase, 
and then dissipate within 3 Myr. This accounting im- 
plies GMC lifetimes of approximately 20to30 Myr. In 
support of the claim that the GMC classification scheme 
constitutes an evolutionary sequence, the authors note 
that among the resolved GMCs in the NANTEN survey, 
Type 3 GMCs are on average more massive, have larger 
turbulent line widths, and have larger radii. However, 
there is significant scatter in the Type 3 GMC sample 
and the mass and size evolution are well within their 
error bars. 

In order to correct for extinction, which might obscure 
Ha emitting H n regions, radio continuum maps at three, 
well-separated frequencies were used to identify obscured 
H II regions via their flat spectral slopes. However, no 
H II regions were identified in the radio continuum data 
that were not present in the Ha maps, leading the au- 
thors to conclude that the Ha data was unaffected by 
obscuration. No similar analysis was performed to esti- 
mate obscuration of young star clusters. No attempt was 
made to correct for the varying sensitivities in the dif- 
ferent radio maps, allowing for the possibility that some 
H II regions were detected at 1.4 Ghz but below the sen- 
sitivity limit at 4.8 and 8.6 Ghz. 

There are several observational biases inherent in the 
GMC classification scheme described above. The first is 



the probable existence of star clusters and H n regions 
located either behind or within giant molecular clouds 
from our viewpoint. High dust extinction along these 
sightlines would mask some young clusters from detec- 
tion in the iBica et all ()1996f ) star cluster sample. This 
could lead to an overestimate of Type 2 GMCs relative 
to Type 3 GMCs. A nother possible bias is the use of 
the IBica et ail (|1996t ) star cluster catalog. Clusters in 
this catalog were targeted for UBV photometry based 
on brightness and association with emission nebulae. It 
is possible that some young clusters were missed in this 
catalog and no attempt is made by Kawamura et al. to 
correct for the completeness of the cluster catalog. This 
would also lead to an overestimate of Type 2 GMCs rel- 
ative to Type 3 GMCs. 

In order to make a quantitative comparison be- 
tween our model s and the evolutionary classification of 
IKawamura et ahl (j2009f ). we employ a few simple pre- 
scriptions to generate synthetic V band, Ha, and ra- 
dio continuum photometry for the clusters and H II re- 
gions in our simulated clouds. First, we calculate the 
U-band luminosity of o ur simulated clusters usin g the 
synthetic photometry of lLeieune fc Schaererl pOOlh . For 
5 M < M* < 15 M^, the photometry i s based on the 
evolution tracks of iSchaerer et all (|199~3l ). For massive 
stars, M* > 15 M Q , the synt hetic photometry is based 
on the high mass-loss models of lMevnet et al.l (I1994D . For 
both sets of synthetic photometry , we assume Z = 0.008. 
Following iParravano et ahl ()2003l ). we approximate vari- 
ations in Ly for these stars by only considering the main 
sequence evolution and taking Ly = (Ly} MS , the mean 
luminosity on the main sequence. Since the stellar evo- 
lutionary tracks give the luminosity at a discrete set of 
masses, we interpolate by fitting a broken power law be- 
tween stellar masses with evolutionary tracks. 

We calc ulate the Ha luminosity o f our model H II re- 
gions via (|McKee fc Williamslll997l) . 

Lhq = 1.04 x lO 37 ^ erg s _1 , (39) 

where 649 is the ionizing luminosity of the central star 
cluster in units of 10 49 photon s _1 . This is larger than 
the empirical relation by a factor of 1.37 to correct for the 
absorption of ionizing radiation by dust grains. Lastly, 
we find the radio continuum lu minosity of our simulated 
H II regions via (|Condonlll99"l . 

L v = 1.6 x 10 23 54 9 (^7^) cr S s_1 Hz_1 ( 40 ) 

We also account for the reduction in flux fro m GMCs that 
would be larger than the beam size used by iDickel et al.l 
(|2005|) by performing a geometric correction and assum- 
ing that all H 11 regions are placed at the center of the 
model beam. 

To assign our simulated GMCs an evolutionary classi- 
fication, we extract the ionizing luminosity and V band 
luminosity of the brightest cluster in the GMC as a func- 
tion of time. Using the ionizing flux, wc calculate the 
expected Ha and radio continuum luminosity via equa- 
tions (|3^)l and (|4TJ|) respectively. For the V band and 
Ha luminosity, we correct for t he foreground extinc tion, 
Ay = 0.25, towards the LMC (jSchlegel et al.lll998ft . ig- 
noring extinction internal to the LMC but external to 
the cloud under consideration. To identify the GMC as 
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Figure 7. Fractio n of GMC lifetime spent as a Type 1, 2, and 3 
GMC as denned by IKaw amura ct al. (2009) (diamonds, top row), 
and average age of GMCs in each classification bin (bottom row). 
The asterisks indicate the median of the GMC type probability 
distribution functions generated using a Monte Carlo analysis de- 
scribed in the text. The error bars encompass the 10th to 90th 
percentile interval of the probability distribution functions. 



either Type 2 or Type 3, we use either an Ha luminos- 
ity cutoff or a radio continuum flux c utoff correspond- 
ing t o the detection limits quoted by (|Kawamura et al.l 
l200flh . For Ha, we say that an H II region is detected if 
Lho > 10 36 erg s _1 . For radio continuum, we say an H n 
region is detected if the radio flux at a distance of 50 kpc 



would be greater than 0.7 mJy bea m 1 for a beam size of 
20" at 4.8 GHz picket et alj|2005l ). Finally, we say that 
a GMC is Type 3 if it meets the criteria just described for 
Ha or radio continuum as well as i f Ly > 1.66 x 10 4 L Q , 
the completeness limit quoted by iBica et all (|1996l ) for 
the young cluster sample. 

Since there are bound to be clusters that are located 
both behind and within clouds along our line of sight, we 
also correct the synthetic V^-band and Ha photometry 
for extinction by the GMC. This is done by assuming 
that star clusters form at random locations within clouds. 
We calculate the optical depth to the star cluster via 
t v = Ku^Rd where K v is the mean dust opacity through 
the cloud, £ M is the depth into the cloud in units of the 
cloud radius and (i is the angle between t he normal to the 
surfa ce of the cloud and the line of sight (|Krumholz et al.l 
I2008T ) . For this purpose, we use a Milky Way extinction 
curve and assume an LMC dust-to-gas ratio whereby a 
column of 1 g cm -2 corresponds to Ay = 107. In the 
visual passbands we are concerned with here, the Milky 
Way and LMC extinction curves arc nearly identical. 

We have run a set of 2000 cloud models, 1000 each 
evolved using two different values of the ambient surface 
density, S rcs = 8 and 16 M pc -2 . All other parameters 
are as in Table [TJ The resulting cloud models encompass 
the entire observe d range in cloud masses reported by 
iFukui et al.l (12001 . 

To directly compare to the observed distribution of 
GMC types we perform simulated observations using a 
Monte Carlo scheme. Since the observations are inher- 
ently weighted by the GMC mass function, we first gen- 
erate cloud masses by drawing from a powerlaw GMC 
mass spectrum with a slope of -1.6, a minimum mass 
of 5 x 10 4 M w and a ma ximum mass of 5 x 10 6 M Q 
(jFukui fc Kawa mura 2010). Once a mass is generated, 
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Figure 8. GMC classification as a function of time for a selection 
of model clouds. 



we find all time steps where model clouds have masses 
within 0.1 dcx of the randomly selected mass. Within 
this sample of time steps, we calculate /i, /2, and fa, 
the fraction of Type 1, Type 2, and Type 3 GMCs, re- 
spectively. At the same time, we calculate ti, t^, and £3 
the average age of clouds in each GMC Type bin. We 
generate 10 4 Monte Carlo realizations, from which we 
construct probability distributions for f±, /2, fz, t\. t-j, 
and t^. 

The results of this comparison are presented in the top 
row of Figure [7] In the figure, the lines connect the me- 
dian of the Monte Carlo probability distributions while 
the error bars encompass the 10th and 90th percentile. 
We are able to reproduce the observe d distribution of 
Type 1 , 2 and 3 GMCs as observed by IKawamura et al.l 
(|2009[) . In particular, using both detection limits, we find 
that the majority of clouds are detected as Type 2 GMCs, 
and relatively fewer clouds are detected as Type 1 and 3 
GMCs. Interestingly, in the bottom panel of the figure, 
we find that, on average, the GMC classification scheme 
does constitute an age sequence in that Type 2 GMCs 
tend to be somewhat older than Type 1 clouds. Type 3 
GMCs in turn tend to be older than Type 2 clouds. On 
the other hand, the spread in cloud ages within each bin 
is well within the error bars, indicating that the GMC 
type classification is not necessarily a strict evolutionary 
sequence: older clouds can be classified as Type 1 and 
younger clouds can be classified as Type 3. 

This can be seen more clearly in Figure [8j where we 
plot the classification of a selection of GMCs as a func- 
tion of time. We see that for some clouds, the classifica- 
tion scheme does represent an evolutionary sequence. In 
these runs, thw cloud starts as a Type 1 GMC, begins 
forming star clusters, evolves into a Type 2 GMC, and 
then forms massive OB association and becomes a Type 
3 GMC. However, we also see that a cloud can quickly 
form a massive OB association and be classified as a Type 
3 GMC early on and only later be classified as a Type 
2 GMC. Alternatively, a cloud may happen to not form 
any massive clusters late in its evolution, causing a mas- 
sive and old cloud to be identified as a Type 2 or 1 GMC 
late in its evolution. Finally, there are clouds which ex- 
hibit no discernible pattern in their histories, more or 
less randomly transitioning between GMC classifications 
throughout their lives. This may explain the presence of 
massi ve (~ 10 6 M^) Type 1 " youn g" GMCs in the sam - 
ples of lKawamura et al.l (|2009D and lHughes et al.l (|2010D . 

6. CAVEATS AND LIMITATIONS 
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6.1. Implications of the Assumption of Homology 

Assuming that clouds evolve homologously is the main 
limitation of our model. We make the assumption of ho- 
mology to significantly simplify the equations governing 
the evolution of the cloud. Given the assumption of ho- 
mology, our equations of motion follow rigorously from 
the local equation of momentum and energy conserva- 
tion. A more complex cloud structure destroys the rela- 
tive simplicity of the model and would require computa- 
tionally expensive hydrodynamical simulations to model 
accretion in detail. 

Homology constrains the cloud to always maintain the 
same shape and degree of central concentration. This is 
equivalent to setting time derivates of the structure con- 
stant, ai, to zero. This might be a problem if changes of 
the moment of inertia of clouds occur primarily by chang- 
ing the shape of the cloud rather than through overall ex- 
pansion or contraction. It might also be a problem if the 
accretion flow does not in reality differentially mix with 
the cloud. If the accretion flow is anisotropic or if it can- 
not fall to the central regions of the cloud before it mixes, 
the cloud may become less centrally condensed and we 
may overestimate the kinetic energy injection since ma- 
terial cannot fall to the bottom of the cloud potential 
well. 

While we cannot resolve the dynamical effect of 
changes in the shape of the cloud, we can resolve changes 
in the size of the cloud. Given the observed Larson 
scaling relations, we expect more massive clouds to be 
larger, implying that clouds must expand as they accrete 
mass. If clouds do form via gravitational instability, they 
must accrete significant mass. We do resolve this behav- 
ior in our models. We caution that our virial analysis 
most readily describes a relaxed system. We may be do- 
ing a poor job of resolving phenomena that occur on a 
timescale comparable to the crossing time. 

6.2. Magnetic fields in the Atomic Envelope 

One key assumption we made in deriving the global en- 
ergy equation was that the atomic envelope contributes 
negligibly to the total magnetic energy associated with 
the cloud. That is, we did not includ e th e magnetic field 
when calculating £ a mb in equation ()C6|) . The motiva- 
tion for this assumption is based on comparisons of ob- 
servations of magnetic fields in dense molecular clumps, 
where it is possible to measure the magnetic field di- 
rectly via the Zeeman effect in lines of OH a nd CN 
(|Troland fc Crutcherll200l iFalgarone et al.ll2008D . to the 
magnetic fields in the atomic ISM, measured via the 
Zeeman effect in the 2 1 cm line of neutral Hydrogen 
(jHeiles fc Trola nd 2005). These studies consistently find 
that the magnetic field strength is significantly elevated 
in the dense molecular gas. However, most of the volume 
of a GMC is occupied by diffuse molecular gas with low 
OH abundance — frustrating efforts to measure magnetic 
fields via observations of OH in emission. Thus, without 
direct observations of magnetic fields in the diffuse molec- 
ular gas, we cannot know whether the magnetic fields 
in the atomic envelope are weak or strong compared to 
the magnetic field strength in the bulk of a GMC. While 
preliminary observations by Troland et al. (2011, private 
communication) indicate that magnetic fields in the dif- 
fuse molecular component are somewhat stronger than in 



the atomic envelope, this question has yet to be settled. 
It is possible that we underestimate the net magnetic 
energy due to the presence of the cloud. 

Another possible problem with our treatment of mag- 
netic fields is that we assume the mass-to-flux ratio re- 
mains constant throughout the evolution of the cloud. 
This is equivalent to assuming a fixed value for rjs- 
This assumption might be invalid if accreted mate- 
rial flows preferentially along magnetic field lines or 
if ambipolar diffusion can act act on timescales com- 
parable to the cloud dyna mical time. While mea- 
surements (ILi et al. 1200611) and theoretical esimatcs 
(|McKee fc Ostrikerl I2007TT of the mass-to-flux ratio of 
molecular clouds find that clouds should be marginally 
supercritical, with mass-to-flux ratios of order unity, the 
time-evolution and cloud-to-cloud variation in mass-to- 
flux ratios are poorly constrained. 



6.3. Validity of Larson's Laws 



It has been suggested 
stancy of GMC su rface 
(see e.g. iKeeell 
1997t iBallesteros 



that the 
densities 



observed con- 
is an artifact 



1989: IVazauez-Semadeni et al.l 
-Paredes fc Mac Lowl~ 120021: 
Ballesteros- Paredesl I2006E Ballesteros-Parcde s et al.l 
20111) . Since the physical structure we assume for the 
clouds in our model implicitly assumes that the Larson 
relations arc valid, this may imply that our models do 
not correspond well to real GMCs. 

The argument that the Larson laws are a product of a 
selection effect usually proceeds as follows. If one looks 
for overdensities of a particular size in a simulation of the 
turbulent ISM, one will find that the selected clouds have 
a wide distribution of masses, implying a large spread 
in cloud-to-cloud surface density. Since observers detect 
clouds using CO as a tracer and at low surface density 
the CO abundance is not high enough to be detectable 
in emission, observers will never find low surface density 
clouds. This implies that at a fixed radius, real clouds 
should have more variation in mass than a naive inter- 
pretation of the CO observations would suggest. 

This argument misses two key aspects of the observed 
properties of GMCs. The first is that it cannot explain 
the lack of clouds at high gas surface densities. For the 
same reason that we should not be able to see diffuse 
clouds, we should very easily be able to see compact, high 
surface density CO clouds. The fact that these clouds 
don't exist implies something important about molecu- 
lar cloud structure. The second argument is that the lack 
of low surface density clouds does not imply that molec- 
ular clouds can exist at all surface densities but merely 
that the molecular clouds dissociate once they become 
optically thin to the ambient ultraviolet radiation field. 
While diffuse ato mic clouds certainly exi st, these clouds 
do not form stars (jKrumholz et al.ll2011| ). 

Since high surface density GMCs are not observed in 
the local universe and low surface density clouds are 
not molecular and thus not GMCs by definition, the ob- 
served lack of variation in GMC surface densities must 
be a real property of the clouds. This has been con- 
firmed with very detailed dust extinction measurements 
of nearby star forming clouds, wh ere an exquisitely tight 
mass-radius relation is observed ([Lombardi et al.ll2010D 
and in extragalactic studies where little variation is seen 
when comparing the mass-radius relation from galaxy 
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to galaxy (jBolatto et al.l [20081) . Taken together, the ev- 
idence seems to imply that the Larson relations are a 
property of the structure of GMCs and are not due to a 
selection effect. 

7. CONCLUSIONS 

In this paper, we have presented semianalytic dynami- 
cal models for the evolution of giant molecular clouds un- 
dergoing both mass accretion and star formation. These 
models are able to capture the evolution of individual 
GMCs from their growth and the onset of massive star 
formation, until their dispersal via an energetic H II re- 
gion or through the combined action of accretion and 
star formation. We are able for the first time to syn- 
thesize galactic populations of GMCs whose properties 
correspond closely to the observed properties of GMCs 
in the Milky Way and nearby external galaxies. We have 
shown that clouds in low surface density environments 
generally disperse within a few crossing times, before 
they can accrete all of the gas in their reservoir. At 
the same time, clouds in high surface density environ- 
ments do accrete all of the gas in their reservoirs and 
tend to be larger and more massive. We have also shown 
that mass accretion can contribute a significant frac- 
tion of the total energy available for turbulent driving. 
Lastly, we generate synthetic cluster observations and 
co mpare against the evo lutionary classification scheme 
of lKawamura et ahl ()2009l ) , finding good agreement when 
we correct for selection effects and systematic biases in- 
herent in the observations. We conclude that, on aver- 
age, the evolutionary classification scheme corresponds 
to an age sequence but is not a good predictor for the 
evolutionary state of isolated clouds. 
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APPENDIX 

A. DERIVATION OF THE VIRIAL THEOREM FOR AN ACCRETING CLOUD 

Here we derive an equation governing the virial balance of a c loud that is simultaneous ly forming stars and accreting 
material. This is a generalization of the analysis of Paper I and Mc Kee fc Zweibel (|1992T ) . We refer the reader to those 
papers for details that are unrelated to the accretion flow and to § 12.11 for a general overview of the model. 

Consider a single molecular cloud contained within an Eulerian volume Kdr with bounding surface SVir- We assume 
that Vvir is sufficiently large to contain the cloud at all times. Material within the virial volume is apportioned into 
three components: virial material, a gaseous reservoir, and material in a photoionized wind. Locally, each component 
satisfies its own continuity equation, 



^ = -V-pv + p (Al) 

-^p = - V • PresVres - Pace (A2) 

= -V ■ PwV w - Ay- (A3) 
These equations are coupled via the source and sink terms, 

P = Pacc+Pcj- (A4) 

We assume that accretion can only transport material onto the cloud and that the wind can only carry material away 
from the cloud, implying p acc > and p e j < 0. 
The local equation of momentum conservation for the virial material is (c.f. equation [A2] in Paper I) 

d 

— (pv) = -V • (LT - T M ) + Pg + F» + Pej(v+ <j) + paeeV rcs (A5) 

where II = P t hl + pvv — tt is the gas pressure tensor, tt is the viscous stress tensor, Tm = [BB — (l/2)_B 2 I]/(47r) is 
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the Maxwell stress tensor, B is the magnetic field, I is the unit tensor, g is the gravitational force per unit mass, and 
F* + /0 C j(v + Vgj) is the local body force due to the interaction between expanding H n regions and the cloud. We 
write the stellar forcing term in this form so as to separate the random component, F», from the spherically symmetric 
component p c j(v + Vgj). 

F* is a function of time and position in the cloud and depends on the precise location and history of massive star 
formation. Due to the nature of supersonic turbulence, we expect turbulent overdensities to be randomly scattered 
throughout the cloud. Thus, we expect F* to be randomly oriented with respect to the position vector, implying 
J v r ■ F*dV = 0. While F* is randomly oriented, p C j(v + v' e -) should on average be purely radial because the 
photoionized wind will be blown out preferentially along the pressure gradient. Neither F* nor the viscous stress 
tensor is included in the momentum equation used in Paper I. A detailed comparison of the two approaches leads us to 
conclude that the formulation used here more properly follows the transport of energy through the turbulent cascade. 

After taking the second time derivative of the cloud moment of inertia, I c \ = J v pr 2 dV, and substituting the 

momentum equation (|A5|I into the resulting expression, we find 
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Upon evaluating the integrals in equation (|A6[) term by term, we obtain the Eulerian Virial Theorem, 

i/ cl = 2(T - To) + B + W - \j t [ (pvr 2 ) ■ dS + ai M cl R cl R cl 
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+ - ai M cl R^ + aiM ej R cl R cl 



s viI 

3 kr 



R c l(M e - 3 v' e : - ^MaccWosc) 



(A7) 



Here T, To, B, and W are respectively the standard kinetic, surface kinetic, magnetic, and gravitational terms (see 
Paper I for precise definitions) and oi = (3 — k p )/(5 — k p ). The final term in equation (|A7[) does not appear in the 
EVT derived in Paper I and is due to the presence of the accretion flow. This has the same form as the wind recoil 
term except for the presence of a dimcnsionless factor 



= J (4 - k p )x 3 - k " { 



fir 
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dx 
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The 



which arises because material is accreted at a velocity that depends on the depth of the cloud potential well, 
dimcnsionless variables x, y, and / are defined explicitly in Appendix [B] 

The magnetic term B retains the form used in Paper I and derived by (jMcKee h Zweibellll992l ) because we assume 
any deformation of the magnetic field in the ambient medium caused by the presence of the cloud is negligible at the 
virial surface, allowing us to approximate that the virial volume is threaded by a constant magnetic field Bo- Here Bo 
is the RMS value of the true magnetic field at the virial surface, which may fluctuate around Bo- Since material in the 
reservoir should also carry currents and thus generate magnetic fields, this parameterization underestimates the total 
magnetic energy. However, we expect the mean density of reservoir material within the virial volume to be an order 
of magnitude smaller than the density of cloud material (see Appendix |5J , so the contribution of the reservoir to the 
magnetic energy is small compared to the contribution due to the cloud. 

B. PROPERTIES OF THE RESERVOIR 

Consider an accreting cloud that is not forming stars and thus not generating a wind. Let M res (r, t) be the mass of 
material in the accretion flow contained within a radius r at time t and let Ar = v TeStSys A.t be the distance that the 
accreting gas falls in a time At. In the same time, a fraction of the material in the accretion flow will be converted 
into cloud material. In the frame comoving with the accretion flow, the change in the mass of reservoir interior to a 
radius r in a time At is, 



M res (r, t) - M ies (r -Ar,t+ At) = -At / Pacc dV. 



(Bl) 

Upon Taylor expanding in At and Ar, dropping the nonlinear terms, and evaluating the integral on the right hand 
side of equation (|B3[) . we find, 



Ar 



dM Tes (r,t) dM Ies (r,t) 



dr 



dt 



-M acc (r,t)At. 



(B2) 



If the accretion flow is in quasi-steady state, the time derivative vanishes and M res (r, t) — M Tes (r). Integrating to 
obtain M Tes (r), we find 
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Figure 9. (a): y (dotted line) and z (dashed line) as a function of x = r / R c \. (b): y(l) = M TeB (R c i)/M c i as a function of £ = M aC ciff/M c i. 
We see M res (i? cl ) > M cl for £ > 5. 
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If we wish to evaluate the above integral and obtain an expression for M res (r), we need to know v ros . sys (r). Expanding 
equation <(4]), we see 
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Equations (|B3[) and (|B4|) constitute a system of integral equations that must be simultaneously solved to obtain a 
solution for the velocity and mass profile of reservoir material within the cloud volume. If we define the functions 
x = r/R c \, y(x) = M rcs (R c ix) / M c \, and / = M C \/(M ICS + M c i), equation (|B4[) reduces to 
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Defining ( = M acc tff/M c \ and = f y(x')/x' 2 dx' , we sec that the problem of determining both 
M rcs (r) reduces to solving the following system of nonlinear ordinary differential equations: 



(B5) 
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(B6) 
(B7) 



These equations can be solved numerically using the shooting method as follows. Since there should be no unmixed 
accreted material at r = 0, we expect y(0) = 0. Since z(x) is defined in terms of an integral, we have no a priori 
knowledge of 2(0). However, we do know that z(l) = and we expect y(l) < £. If we impose y(l) = cQ where c is a 
constant of order unity, we can integrate the above system and obtain a trial solution for y{x) and z{x). The constant 
c can then be varied until a solution is obtained with y(0) = 0. There is a unique solution for each value of £ and thus 
a new solution must be obtained if £ varies. 

A key assumption in this analysis was that the accretion rate is approximately constant over the cloud free-fall 
timcscale. This is equivalent to the assumption that £ < 1. This is a reasonable assumption, which we can see by 
making an analogy to the case of a protostellar core. Since we expect M c \ « M c i/tg tI where tg tI is the free-fall time for 
all of the gas in the reservoir UStahler et al.|[l980t iMcKee fc Tar][200l and since the reservoir should be significantly 
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more extended than the cloud, we expect the mean density of the reservoir to be much lower than the mean density 
in the cloud. This implies ts,r 3> iff el- If this condition does hold, then the condition ( < 1 is automatically satisfied. 

In Figured we present a numerical solution for y(x) and z(x) obtained for £ = 1.0. We sec for this case that 
M res (i? c i) w 0.2M c i, showing that even for substantial accretion rates, the gas within V c \ is primarily composed of 
cloud material. In FigurelHt) we show y(l) = M ros (P c i) / M c \ as a function of £. Reservoir material becomes the primary 
component of the volume occupied by the cloud for £ > 5.0. For this reason, we reject as unphysical any portion of 
cloud history with ( > 5. In practice, we find ( < 1.0 for the lifetime of all clouds we simulate. 

C. DERIVATION OF THE EQUATION OF ENERGY CONSERVATION FOR AN ACCRETING CLOUD 

Here we derive the global equation of energy conservation for a cloud undergoing accretion and star formation. We 
begin by writing the equation of momentum conservation (equation |A5| ) in Lagrangian form, 

= -VP - P V4> + V • 7T + + PcjV^j + F„ + Pacc(v rcs - v). (CI) 

Upon contracting the above equation with v, we obtain the nonthermal energy evolution equation, 



d ( 1 B 2 — B 2 \ ( \ 

— f -pv 2 + p0 c i + — 2- J + V • pv ( -v 2 + cj) cl J + V • S P = -v ■ VP + p^- + pv ■ g ros + V • (tt ■ v) 



(C2) 

7r : Vv + p f -v 2 + <j) cl J + v • F* + p cj v • + p acc (v ■ v rcs - v 2 ^ 

where Sp is the Poynting vector. Combining the first law of thermodynamics with the continuity equation yields the 
evolution equation for the thermal energy of the cloud, 

^(pe)+V ■ pv(e+^) =pe + v VP + T-A (C3) 

where e is the internal energy per unit mass and T and A are respectively the rates of energy gain and loss per unit 
volume. Here we've assumed that accreted material has the same thermal energy density as cloud material, which will 
be true if the radiative timcscalcs are short compared to the mechanical timescales, as in molecular cloud conditions. 
Summing equation (|C2j) and (|C3j) , we obtain the evolution equation for the total energy of the system 
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p [ -^v + e + -4>d j + -p<pci + pv • g res + pejv • v oj + p acc (v • v rcs —v ) + V • (tt ■ v) + v • F* - A. 



Here we've assumed that r = it : Vv, the scalar rate of viscous dissipation. This is equivalent to the statement that 
turbulent kinetic energy is converted into heat at viscous scales whereupon the energy is quickly radiated away. While 
there may be other heating mechanisms, including cosmic rays and protostellar radiation, we neglect these sources by 
noting that any local heating should be offset in a time much shorter than the dynamical time. 
If we define the total energy due to the presence of the cloud, 

£ * = I P (l v " + e + U*) dV + I ^o^ dV (C5) 



JV cl \ 2 2 / ./Kir 8lT 

and the total energy in the ambient medium, 

£amb = ^ ^ p(^v 2 + e+^ c ijdV, (C6) 

we can integrate equation (|C4j) over the virial volume and obtain an evolution equation for £ c \. Portions of this 
calculation arc performed explicitly in Paper I and we do not reproduce those results here. The new terms stem 
from our inclusion of accretion as well as the inclusion of F», the random component of the stellar forcing term. The 
integrals over the new terms are evaluated below. 
The first new term is due to the gravitational influence of the reservoir, 

P c i GM c] 

pv • g rcs dV = -x-5 ^ (C7) 

where 
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X = (3-fc p ) / x 1 - k ^y{x)dx. (C8) 
Jo 

This is the gravitational work done on the cloud by the reservoir material as the cloud expands and contracts. 

The random stellar forcing term can be evaluated by noting that r and F* arc assumed to be uncorrelated, so 
J v v • F*dV = J v Vt ur b • F*dV. This integral depends on the degree of correlation between two randomly oriented 
vectors, v tur b and F„. Since, at a fixed time, the direction of expansion of an H II region is the same as the direction of 
momentum injection, these vectors should indeed be highly correlated. This term then corresponds to the net injection 
of energy by H II regions, 

Gel = ! Vturb • F*dV. (C9) 

"'Kir 

For the terms proportional to p acc , we have 



and 



PaccV • VresdV = - [ 7—^ ) £,RclM acc V csc + / Pace Vturb ' V rcS dV (CIO) 
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p acc v 2 dV = aiM acc R 2 cl + 3M acc CT c 2 ! (Cll) 



Kir 

where we've used our assumption that v re s,rand is randomly oriented with respect to r. 

We are left with one last integral on the right hand side of equation (|C10p that depends on the correlation between 
two vector fields, v tur b and v ros . Suppose v rcs and v tur b are perfectly correlated. In that case, the integral we are 
concerned with can be readily evaluated in three limits, |v rcs | = |v tm b|, |v ros | 3> |v tm b|, and |v res | -C |v tU rb|. In the 
first limit |v res | = |v tU rb| the transfer of material to the cloud is merely an act of relabeling, so v tU rb • v res = v^. Now 
consider the limit |v res | 3> |v tU rb|- We approximate that a parcel of reservoir material mixes with the cloud once it has 
swept up a mass of cloud material equal to its own mass. Since the interaction between cloud material and reservoir 
material must be inelastic if they are to mix, the velocity of cloud material must be driven by the act of mixing to 
v tmb = v rcs /2. Conversely, if |v res | <C |v t urb|: the reservoir material is driven to a velocity v rcs = v tU rb/2- We obtain 
the correct answer in all three limits if we assume 

PaccVturb • V Ies dV = i jf pacc^L + V? urb )dV. (C12) 

This is an upper bound on the value of the integral in the limit of perfect correlation between v tur b and v res . If v tur b 
and v res are uncorrelated then the integrand should on average be zero. Thus we have, 

Pace Vturb ' V res ^ < ^ J Pacci^res + v tuib) dV - ( C13 ) 

To evaluate this integral, we linearly interpolate between the upper limit and lower limit, 

If 3 • 3 • 
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where ip is an interpolation parameter that ranges between zero and unity and 

7 = U/( / T7T 1 7T--^7-V (C15) 
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Summing the individual terms derived above yields the global form of the equation of energy conservation, 

(C16) 
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where £ c i = J v AdV. 
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